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The general method for treating non-Gaussian wave functionals in the Hamiltonian formulation 
of a quantum field theory, which was previously proposed and developed for Yang-Mills theory in 
Coulomb gauge, is generalized to full QCD. For this purpose the quark part of the QCD vacuum 
wave functional is expressed in the basis of coherent fermion States, which are defined in term of 
Grassmann variables. Our variational ansatz for the QCD vacuum wave functional is assumed to be 
given by exponentials of polynomials in the occurring helds and, furthermore, contains an explicit 
coupling of the quarks to the gluons. Exploiting Dyson-Schwinger equation techniques, we express 
the various n-point functions, which are required for the expectation valnes of observables like the 
Hamiltonian, in terms of the variational kerneis of our trial ansatz. Finally the equations of motion 
for these variational kerneis are derived by minimizing the energy density. 
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I. INTRODUCTION 

One of the major challenges of theoretical particle physics is the understanding of the low-energy sector of Quantum 
Chromodynamics (QCD). Despite many years of intensive research, a thorough and unified picture of the low-energy 
phenomena of strong interactions, i.e., confinement and spontaneous breaking of chiral symmetry, is still lacking. 
Much insight has been gained by means of lattice Monte-Carlo calculations, in particular on the gluon sector of QCD. 
However, despite much progress, the treatment of dynamical chiral quarks is still a challenge for the lattice approach, 
which furthermore struggles to describe the phase diagram of QCD at finite baryon density due to the notorious 
sign problem. In addition and in general, a thorough understanding of physical phenomena cannot be achieved by 
numerical lattice simulations alone but analytic methods, albeit approximate ones, are needed as well. 

Over the last decade substantial efforts have been undertaken to develop non-perturbative continuum approaches 
to QCD. These approaches are based on functional methods and can be roughly divided into three classes: i) Dyson- 
Schwinger equations (DSEs) in Landau [1-3] and Coulomb gauges [4-9], ii) functional renormalization group (FRG) 
flow equations [10-12] and iii) the variational approach to Hamiltonian QCD [13, 14]. These three approaches are 
intimately related. The first two approaches are based on the functional integral formulation of QCD in either Landau 
or Coulomb gauge, while the variational approach has been mainly applied to the Hamiltonian formulation in Coulomb 
gauge but has been recently also extended to the effective action of the functional integral formulation in Landau 
gauge [15]. The equations of motion of the first two approaches are in faet very similar and in a certain approximation 
(replacing the renormalization group scale k in the loop integrals by its infrared value k = 0) the FRG flow equation 
becomes a Dyson-Schwinger equation. The FRG flow equations have been also applied to the Hamiltonian approach 
in Coulomb gauge and results similar to those in the variational treatment were found [16, 17]. Furthermore, the 
DSF techniques can be very advantageously exploited to carry out the variational approach with non-Gaussian wave 
functionals describing interacting quantum helds [18]. In the present paper we use the Dyson-Schwinger equations to 
develop a variational approach to the Hamiltonian formulation of QGD. 

Previous variational studies within the Hamiltonian approach have focused on the Yang-Mills sector and used 
Gaussian-type ansåtze for the vacuum wave functional. This has provided a decent description of both the infrared 
(IR) and ultraviolet (UV) sector in rough agreement with the existing lattice data. In particular, a linearly rising 
static quark potential [19, 20] and a perimeter law [21] for the ’t Hooft loop [22] were found, which are both features 
of the conhned phase. Also the connection to the dual Meissner effeet (an appealing picture of conhnement) was 
established [23]. More recently, the deconhnement phase transition was studied at hnite temperatures [24, 25] and 
the effective potential of the Polyakov loop was calculated [26, 27]. The obtained critical temperature is in reasonable 
agreement with the lattice data. Furthermore, the order of the phase transition was correetly reproduced [27]. In 
the zero-temperature calculation the obtained static gluon propagator agrees with the lattice data in the IR and UV 
but misses some strength in the mid-momentum regime. This missing strength can be attributed to the absence of 
non-Gaussian terms in the vacuum wave functional [18]. 

Generally, Gaussian wave functionals describe quantum held theories in the independent (quasi-)particle approx¬ 
imation, while truly interacting quantum helds possess non-Gaussian vacuum wave functionals. In Ref. [18] a general 
method for treating non-Gaussian trial wave functionals in a quantum held theory was proposed. This method relies 
on Dyson-Schwinger-type of equations to express the various n-point functions of the quantum held in terms of the 
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variational kerneis contained in the exponent of the ansatz for the vacuum wave functional. So far this method 
has been formulated for Bose helds only and was applied to the Yang-Mills sector of QCD using a wave functional 
which includes, besides the usual quadratic term of the Gaussian, also cubic and quartic terms of the gauge held. In 
principle, the approach put forward in Ref. [18] is general enough to deal with any interacting quantum held theory. 
In the present paper we extend this approach to full QCD. The central point will be the treatment of fermion helds 
interacting with Bose (gauge) helds. To exploit the Dyson-Schwinger equation techniques, the second quantization of 
the fermion sector of the theory has to be formulated in terms of Grassmann variables. For this purpose we express 
the quark part of the QCD vacuum wave functional in terms of coherent fermion States. We will formulate the present 
approach to Hamiltonian QCD for general wave functionals but work out the Dyson-Schwinger-type of equations only 
for those wave functionals whose exponent is bilinear in the quark held. 

The QCD vacuum wave functional is chosen as the exponential of some polynomial functional of the quark and 
gluon helds. The coefhcient functions of the various polynomial terms are treated as variational kerneis. By means of 
Dyson-Schwinger equation techniques we express the various n-point functions, needed for the vacuum expectation 
value of observables like the Hamiltonian, in terms of these variational kerneis, which in this context hgure as bare 
vertices. The resulting equations are different from the usual DSEs, which relate the various full (dressed) propagators 
and vertices to the bare (inverse) propagators and vertices occurring in the classical action, and are termed canonical 
recursive DSEs (CRDSEs) in the following. By means of the CRDSEs we express the vacuum expectation value of 
the QCD-Hamiltonian as functional of the variational kerneis of our vacuum wave functional. Minimization of the 
energy density with respect to these variational kerneis results then in a set of equations of motion (referred to as 
‘gap equations’), which have to be solved together with the CRDSEs. 

The organization of the rest of the paper is as follows: In Sec. II we briefly summarize the basic ingredients of the 
formulation of the second quantization in terms of Grassmann variables and present the quark part of the QCD wave 
functional in the basis of coherent fermion States. In Sec. III we derive the general form of the CRDSEs for the static 
Green functions of QCD, assuming a QCD vacuum wave functional which in particular contains the coupling between 
quarks and gluons. In Sec. IV, after introducing the QCD Hamilton operator in Coulomb gauge, we calculate the 
energy density in the vacuum state. The variational principle is carried out in Sect. V, where we derive the equations 
of motion for the variational kerneis of our QCD vacuum wave functional. A short summary and our conclusions are 
given in Sect. VI. Here we also briefly discuss further applications of the approach developed in this work. 


II. COHERENT STATE DESCRIPTION OF THE FERMIONIC FOCK SPACE 

As demonstrated in Ref. [18] for Yang-Mills theory, the use of non-Gaussian wave functionals in the Hamilto¬ 
nian approach can be conveniently accomplished by exploiting Dyson-Schwinger equation techniques known from the 
Lagrangian (functional integral) formulation of quantum held theory. This refers, in particular, to wave functionals 
describing interacting helds. To exploit the DSE techniques in the Hamiltonian formulation of QCD it is necessary to 
represent the quark operators and wave functionals in terms of anti-commuting Grassmann helds. In this representa- 
tion the matrix elements between Fock-space States are then given by functional integrals over Grassmann helds. The 
formulation of the second quantization in terms of Grassmann variables becomes particularly efhcient when coherent 
fermion States are used. Below we will briehy review the basic ingredients of the coherent fermion state representation 
of Fock Space and apply it to Dirac fermions. 


A. Coherent Fermion States and Grassmann Variables 

Consider a Fermi system described in second quantization in terms of creation and annihilation operators 5^,, bk, 
satisfying the usual anti-commutation relations 


{bk,bi}=0={blbj}, {bk,bl}=Ski, 


where the subscripts k, I, ... denote a complete set of single-particle States. Let |0) be the Fock vacuum, i.e. 


fcfcjO) =0. 


The coherent fermion States |C) are dehned as eigenstates of the annihilation operators [28] 


&feic) = aic), 


(1) 
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where the Cfe are anti-commuting (Grassmann) variables {CfciO} =0. The corresponding bra-vectors are defined as 
left eigenstates of the creation operators 

(ci&I = (cia. (2) 

where the operation denotes the involution. We will keep the same symbol ‘=i=’ also for the usual complex conjugation 
of ordinary complex numbers. 

The coherent fermion States defined by Eqs. (1) and (2) can be expressed in Fock space as 

IC) =exp(C -6^)10), (Cl = {0|exp(6-C*). (3) 

Here we have skipped the indices and used the short-hand notation 

r,*-C = Y.vtCk- 

k 

From the representation Eq. (3) one easily finds for the scalar product of two coherent States 

(»?IC) = 

The coherent fermion States form an over-complete basis of the Fock space, in which the unit operator has the 
representation 

l = ydC*dCe-‘;*<|C)(CI, With dCMC^n^CfcdCfe. (4) 

k 

An arbitrary state |<?) of the Fock space can be expressed in the basis of coherent fermion States |C) by taking the 
scalar product 


<^(0 = (Cl<?), (5) 

which can be interpreted as the “coordinate representation” of fermion states with the Grassmann variables interpreted 
as classical fermion coordinates. From the representation Eq. (3) it is clear that the ^(C*) are functionals of the Q., 
which can be Taylor expanded 


'*■(';•) = F S - --CL ’ 

n ki...kn 


( 6 ) 


where the ^ki...kn are complex numbers. The representation of bra vectors is obtained by taking the adjoint of Eq. (5)^ 

(<?|C) = (<?(C*))* = <?*(C) (7) 

and is obviously a function of the variables Cfc- From Eq. (6) we find 


n k\...kn 


The Fock-space states are given by functions of the creation operators /(&^) acting on the vacuum state 

\<P)=f{hk\Q). 

From Eq. (2) we find then the corresponding coherent state representation 

HC) = f{C), 

where we have used (C|0) = 1, which follows immediately from Eq. (3). The scalar product between two Fock states 
\<P) and l'F) is easily obtained in the basis of coherent states by inserting the completeness relation Eq. (4) 

{^\<P) = J dC* dC e-^*< ^*(0 (8) 


^ The adjoint means the involution for the Grassmann variables and complex conjugation for ordinary complex numbers. 
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where we have also used the definitions (5) and (7). Using Eq. (2) and 

m = ^ (Cl 

we find for the action of an operator on a Fock State in the coherent State basis 

{C\oib,b^)\<P) = o{^,C)<i>{C) 

and similarly for matrix elements between States of Fock space 

(^^10(6, b^)\<P) = J dC* dC e-^*< W*{C) 0( afv, C*) (9) 

In the next subsection the coherent fermion state representation of the Fermionic Fock space given above is extended 
to Dirac Fermions. 


B. Coherent-State Representation of Dirac Fermions 

The description of Fermi systems in terms of Grassmann variables outlined above can be immediately applied to 
Dirac fermions once the fermion held is expressed in terms of creation and annihilation operators. 

In this paper we use a compact notation in which a single digit 1,2,... represents all indices of the time-independent 
held. For example, for the quark held we have 


^(l)=<(xi), 


where si denotes the Dirac spinor index, while ki stands for the colour (and possibly flavour) index.^ A repeated 
index implies integration of spatial coordinates and summation over the discrete indices (spinor, colour, flavour). 
Furthermore, we define the Kronecker symbol in the numerical indices to contain besides the usual Kronecker symbol 
for the discrete indices also the S function for the continuous coordinates, e.g. 

(I(1,2)=<I(xi-X2)5,,,,5'=^'=^... 


The usual anticommutation relation for the Dirac held reads 

(2)} =<5(1,2). 


( 10 ) 


Let 


ho(l, 2) = (-ia • V + /3m)(5(xi - Xa) (11) 

denote the Dirac Hamiltonian of free quarks with a bare mass m. This Hamiltonian possesses (an equal number of) 
positive and negative energy eigenstates. We can expand the quark held ■i/'(l) in terms of these eigenstates. Let 'ip±{l) 
denote the part formed from the positive/negative energy modes. Obviously, we have 

^(1) = V’+(i) +V'-(i)- 


For the subsequent considerations it will be convenient to introduce orthogonal projectors A± onto the positive and 
negative energy part of the Dirac held, i.e. 


V'±(1) = A±(1,2)V^(2), 


( 12 ) 


satisfying 


A+(1,2)+A_(1,2) = (5(1,2), 


A+(l, 2)A_(2,3) = 0, A±(l, 2) A±(2,3) = A±(l, 3). 


(13) 
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In the present paper the quark flavour will be irrelevant, but the subsequent considerations do not change when flavour is included. 
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These projectors can be expressed by spectral sums over the positive and negative, respectively, energy modes of the 
Hermitian Dirac Hamiltonian ho{l,2), and satisfy A]j_( 1,2) = A±(l,2). With this relation we find for the adjoint 
fermion operator from Eq. (12) 


4(l)=^t(2)A±(2,l). 


(14) 


From Eqs. (12) and (14) it follows with (13) 

{V'±(1),4(2)} = A±(1,2), {^±(1),V'^(2)} = 0. 

In the bare (free Dirac) vacuum state |0) all negative energy modes are filled, while the positive energy modes are 
empty, implying 


V^+(1)|0)=0 = V'1(1)|0). (15) 

In analogy to Eq. (1) we define coherent fermion States |^) = |^+,5_) by 

Mim = c+(i)io, i’Um) = tm), (le) 

where are Grassmann helds. Since V’±(l) = A±(l, 2)'i/)±(2), the Grassmann helds also satisfy 

e±(l)=A±(l,2)e±(2) 

and the positive- and negative-energy component helds can be assembled into a single Grassmann-valued spinor held 

C(l) = e+(l)+e-(l), e±(l)=A±(l,2)e(2) (17) 

satisfying 


= A±(l,2), ^i^=A±(2,l). 

sum ^ ' <5d(2) 

The coherent fermion States dehned by Eq. (16) have the Fock-space representation [cf. Eq. (3)] 

\0 = exp[C-H(l)V’+(l) +?-(l)V'-(l)] |0) 


(18) 


from which follows 

Using Eqs. (16) and (19) the action of Dirac held operators on fermion States \<P) is expressed in the basis of coherent 
States as 

V <^^+(1)/ (20) 

mHm) = + 

where (^|4>) = is the coherent-state representation of the quark vacuum wave functional \<P), which can be 

interpreted as the “coordinate representation” of the latter. Notice that from the dehnition Eq. (16) of the coherent 
fermion States follows that the quark vacuum wave functional depends only on and 

In analogy to Eq. (9) the matrix element of an operator O between fermionic Fock States |^i) and \<p 2 ) is given by 




I V£} e-^ <!>* , 5^] O ^+ ^] ^>2 , ^_] 


where we have introduced the quantity 


( 21 ) 


M = ?|(i)?+(i)-d(i)C-(i), 


(22) 
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arising from the integration measure of the Grassmann fields [cf. Eq. (8)]. Note that fx, being bilinear in the Grassmann 
variables, commutes with any element of the Grassmann algebra. This quantity can be rewritten with the help of 
Eqs. (13) and (17) as 

Pi = C+(l) [A+(l, 2) - A_(l, 2)]C(2) = ^t(i) 2 ) ^(2), (23) 

where the quantity 

Q-\l, 2) = A+(l, 2) - A_(l, 2) = go(l, 2) (24) 


represents the free static quark propagator 

Qo(l,2) = (0|[V^(l),^t(2)]|o), (25) 

with |0) being the Eock vacuum Eq. (15) of the free quarks, as we will discuss in more detail at the end of Sec. HIG. 


C. Coordinate Representation of QCD Wave Functionals 

With the coherent States |^) [Eq. (16)] of the Dirac fermions at hånd we are now in a position to express QGD 
wave functionals in the coordinate representation. The coordinates of the gauge held are its spatial components 
A“(x). In analogy to the fermion held we will collect all indices of the gluon held in a single digit A{1) = A“^(xi), 
where ai denotes the colour index of the adjoint representation and ti is a spatial Lorentz index. In the coordinate 
representation of the wave functional of the Yang-Mills sector tf'[A] = (Al'F) we have for the operators of the canonical 
variables 

{A\Å{im=Ail)nA], (Aln(l)ltf') = 1^. 

We will work here in Coulomb gauge diAf = 0, where only the transverse components of the gauge held are left. In 
our compact notation we have 


where 


M(I) 

6A{2) 




(26) 


t{l,2) = 



^11 


S{xi 


X2) 


is the transverse projector. 

The vacuum state of QCD can be written in the coordinate representation (i.e., coherent-state representation for 
the fermions) as 


=: exp|-i 5^[A] - A]|, (27) 

where Sa dehnes the vacuum wave functional of pure Yang-Mills theory, while S / dehnes the wave functional of the 
fermions interacting with the gluons. On general grounds Sf contains only even powers of Grassmann variables so 
that this quantity, as well as the vacuum wave functional !f'[A,^)i_,^_], commutes with any Grassmann held. 

A comment is here in order concerning the representation Eq. (27) of the vacuum wave functional for the gluon 
and quark helds. Eor the bosonic gluon held we use here the usual “coordinate” representation regarding its spatial 
components as the (classical) coordinates of the theory. As already discussed in Sec. IIA, the classical analogues of 
the fermion helds are the Grassmann variables and the “coordinate representation” of the fermion wave functional is 
the coherent (-fermion) state representation, Eq. (5). We could also use bosonic coherent States for the gluons but 
this is not necessary and we will not use it since the usual coordinate representation is in this case quite convenient. 

Eor sake of illustration we quote the expressions for the perturbative QCD vacuum state [29]. The perturbative 
Yang-Mills vacuum is obtained by choosing in Eq. (27) the quadratic gluonic “action” 

Sa = I ^A“(p)[p|A“(-p). 


(28) 
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On the other hånd, the ground state |0) [Eq. (15)] of the free (perturbative) quarks reads in the coherent State basis 
= (^|0) = 1 (i.e., = 0). This is because the kinematics of the free fermions is already encoded 

in the integration measure fx of the Grassmann helds, see Eqs. (21)-(23). 

Eventually we are interested in the Hamiltonian formulation of QCD in Coulomb gauge. The coordinate represent- 
ation of gluonic States and matrix elements has been presented in Ref. [18]. Writing the QCD vacuum wave functional 
in the form (27) and using the representation (21) of the fermionic matrix elements, expectation valnes in gauge-hxed 
QCD are given by^ 

(0[V>,V’^A^])= j VS}ViVAJAe-^^ + + (29) 

where 

Ja = Det Ga^ (30) 

is the Faddeev-Popov determinant. In Coulomb gauge the Faddeev-Popov operator reads 

= {-5<^^dl-gr^A‘l{^)d-)5{^-y). (31) 

Here g is the bare coupling constant and are the structure constants of the 5u{N) algebra. The functional 
integration in Eq. (29) runs over transverse gauge held conhgurations and is, in principle, restricted to the hrst 
Gribov region or, more precisely, to the fundamental modular region. 

Once the functional derivatives in Eq. (29) are taken, the vacuum expectation value of an operator boils down to 
a functional average of the form 

= J V^^V^VAJa f[AJJ% (32) 

where / is a functional of the helds only (i.e., / contains no functional derivatives). 

For the variational approach to QCD to be developed later we need the vacuum expectation values of products of 
gluon operators A, H, and, in particular, quark held operators ip, Furthermore, to exploit DSEs techniques we 
have to express these expectation values by n-point functions of the Grassmann helds This can be achieved by 

means of Eq. (20). To illustrate how this is accomplished we consider temporarily the quark sector only and omit the 
integration over the gauge held. Using Eqs. (21) and (27) we obtain for the quark bilinear 

(^(l)V't(2))= J v£,^Vi e-'^-^^(^e-(l) + ^^)(?|(2) + ^^) e-^U 

Since S*j is independent of and it is convenient to perform integrations by parts with respect to ^+(1) and 
^_(2).^ The integration by parts with respect to $+(1) yields 

(^(1)V'^(2))= J + e-^f, 

where we have used the hrst equation of Eq. (17). In the same way an integration by parts with respect to C-(2) 
yields 

J e-«?-'^(^^(I)^|(2) + ^^+C(l)^^)e-«/ 

= j (33) 

^ With a slight abuse of notation we will use the same symbol for the gauge held operator and for the field variable to be Integrated over. 

It should be always clear from the context which quantity is meant. 

^ Recall that the formula for integration by parts for Grassmann variables reads 

J dri J dri 
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where again Eqs. (17) and (18) and the normalization of the functional integral were used. Finally, we can generalize 
Eq. (33) to the vacuum expectation valne of full QCD (i.e., taking also the functional average over the gluon sector) 
and include also a functional f[A] of the gauge held. One obtains then 

m) ^t(2) f[A]) = (^(1) ^t(2) f[A]) + A_ (1, 2){f[A]). 

Using the anticommutation relation Eq. (10), from the last relation follows 

(V^t(i) ^(2) f[A]) = (et(l) e(2) f[A]) + A+(2, l){f[A]). 

The expectation value of four fermion operators can be derived along the same lines, resulting in 

(V^t(i) ^(2) ^t(3) ^(4) fiA]) = (et(l) ^(2) ^t(3) C(4) f[A])+ 

+ (^{1) m f[A]) A+(4,3) + (e(2) ^^3) f[A]) A+(4,1) 

+ (^^3) C(4) flA]) A+(2,1) + (e^l) e(4) flA]) A_(2,3) 

+ [A+(2,1) A+(4,3) + A+(4,1) A_(2,3)] {f[A]). 

III. CANONICAL RECURSIVE DYSON-SCHWINGER EQUATIONS OF QCD 

As already discussed in the introduction, in the Hamiltonian approach to a quantum held theory the use of non- 
Gaussian wave functionals is most conveniently accomplished by exploiting DSE techniques known from the functional 
integral approach. In the Hamiltonian approach, the scalar products or matrix elements between the wave functionals 
are given by functional integrals [see Eq. (29)], which have an analogous structure to those of vacuum transition 
amplitudes in the functional integral approach, except for a different action and that the helds to be Integrated over 
live in the spatial subspace only. In general, DSEs relate various propagators (or proper n-point functions) to each 
other through the vertices of the action. In the Hamiltonian approach the “action” is dehned by the ansatz for the 
vacuum wave functional [see Eq. (27)] and the corresponding generalized DSEs, the CRDSEs, are used to express the 
n-point functions, occurring e.g. in the vacuum expectation value of the Hamiltonian or other observables, in terms 
of the variational kerneis occurring in the exponent of the wave functional. 


(34) 

(35) 

(36) 


A. General Form of the CRDSEs 


In the following we derive the CRDSEs for the Hamiltonian approach to QCD. The general structure of these 
equations does not depend on the specific ansatz for the vacuum wave functional. Therefore we will leave the “action” 
^Sa -I- 5/ in Eq. (27) for the moment arbitrary. 

In the preceding section we have seen that any vacuum expectation value of an operator involving the quark helds 
'lp and , and the gluon coordinate and momentum operators A and H, can be reduced to an expectation value of a 
functional of the Grassmann helds and the held “coordinate” A, see Eq. (32). For such expectation values we 
can write down Dyson-Schwinger-type of equations in the standard fashion by integrating a total derivative. 

Taking derivatives with respect to the fermion helds yields 


(,5^(1) 


6S 

TO 



Sf[A,i,e] \ 

sm r 

^e^i) /■ 


where / is any functional of the helds, and 


S = SA + Sf + S} + 


(37a) 

(37b) 


(38) 


is the total “action” dehned by the vacuum wave functional Eq. (27) and with /x, dehned in Eq. (22), arising from the 
integration measure over the fermion helds. 

Similarly, the derivative with respect to the gauge held A leads to 



sf[A,^,e] 

SA{1) 


+ ro(3,2;I)(G^(2,3)/[A,e,C^]), 


(39) 
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where the last term arises from the derivative of the Faddeev-Popov determinant, Eqs. (30) and (31), which gives rise 
to the bare ghost-gluon vertex 


ro(2,3;l) 


5A{1) 


(40) 


Equations (37) and (39) are the basic CRDSEs of the Hamiltonian approach to QCD. By choosing for f[A, succes- 
sively higher powers of the helds one generates from these equations inhnite towers of CRDSEs. The three different 
towers of equations obtained from Eqs. (37) and (39) are equivalent to each other as long as no approximations are 
introduced. 


B. Ansatz for the Vacuum State 

The CRDSEs derived in the previous subsection are quite general, and their structure does not depend on the 
details of the specihc ansatz for the vacuum wave functional Eq. (27). In order to proceed further we have to specify 
the form of the vacuum wave functional. Since we are mainly interested here in the quark sector, we will use a 
Gaussian wave functional for the Yang-Mills sector 

Sa = uj{1,2)A{1)A{2); (41) 

the generalization to non-Gaussian functionals is straightforwardly accomplished by using the GRDSE approach 
developed in Ref. [18] for the Yang-Mills sector. 

The fermionic vacuum state is chosen in the form 

Sf = e+a) KAil, 2) e-(2) = e^(l) A+(l, 1') Ka{1', 2') A_(2', 2) ^(2) (42) 

where the kernel Ka is supposed to contain also the gauge held, so that Sf includes also the interaction of the quarks 
with spatial gluons. We choose this kernel in the form 

KAil, 2) = Koil, 2) + Kil, 2; 3) A(3), (43) 

where Kq and K are the variational kerneis with respect to which we will later minimize the energy. The ansatz 
Eq. (43) can be considered as the leading order expansion of Ka in powers of the gauge held. Furthermore this ansatz 
specihes the quark wave functional as Slater determinant, which guarantees the validity of Wick’s theorem in the 
quark sector. For later use we also quote the complex conjugate fermionic functional 

S} = ^+(1) A_(l, 1') K\il', 2') A+(2', 2) ^(2), 

where 

A^(l, 2) = [KAi2, 1)]* = K*i2, 1) + K*i2, 1; 3) A(3). 

With the ansatz specified by Eqs. (42) and (43) the fermionic part of the action Eq. (38), SfSj -\-can be rewritten 
as 

Sf + S}+fi = C(l) [Qo'(r 2) + 7(1, 2) + fo(l, 2; 3) A(3)]C2), (44) 

where Qo is dehned by Eq. (24) and we have introduced the biquark kernel 

7(1,2) = A+(l, T) Koil', 2') A_(2', 2) + A_(l, T) K^il', 2') A+(2', 2) (45) 

and the bare quark-gluon vertex 

fo(l, 2; 3) = A+(l, T) Kil', 2'; 3) A_(2', 2) + A_(l, T) K^l', 2'; 3) A+(2', 2). (46) 

With the explicit form of the vacuum wave functional given by Eqs. (41) and (44) the CRDSEs (37) and (39) take 
the following form 

[Qo-i(l, 2) + 7(1, 2)] (^(2) /) + fo(l, 2; 3)(C(2) A(3) /) = (47a) 

[Q-\2, 1) + 7(2, 1)] (C(2) /) + fo(2,1;3)(C(2) A(3) /) = (47b) 

20.(1, 2)(A(2) /) - fo(3, 2; l)(e(2) C(3) /) - ro(l; 3, 2)(G^(2, 3) /) = 


(47c) 
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Figure 1. Diagrammatic representation of the CRDSE (52) for the quark propagator. Continuous and wavy lines represent, 
respectively, quark and gluon propagators, where thin and thick lines stand, respectively, for bare and fully dressed propagators. 
Empty boxes represent variational kerneis while full dots stand for one-particle irreducible vertex functions. A dictionary of 
our diagrammatic conventions is given in the appendix. 


Choosing the functional / appropriately, these equations allow us to express the various n-point functions of the 
helds A, in terms of the (variational) kerneis w [Eq. (41)], Kg and K [Eqs. (42) and (43)] of the vacuum wave 
functional. Later on we will use these equations to express the expectation value of the QCD Hamiltonian in terms 
of these variational kerneis. 


C. Quark Propagator CRDSE 

Since our quark wave functional is the exponent of a quadratic form in the fermion helds [see Eq. (42)] we can 
express all fermionic vacuum expectation values of these helds in terms of the fermionic two-point function, which is 
a manifestation of Wick’s theorem. The CRDSE for the full quark propagator 

Q(l,2):=(e(l)et(2)) (48) 

can be obtained by putting / = in Eq. (47a). This yields 

[go-i(l, 3) + 7 ( 1 ,3)] g(3, 2 ) + fo(l, 3; 4)(^(3)^t(2)A(4)) = 5(1, 2 ). (49) 

To resolve the occurring three-point function we introduce the full quark-gluon vertex® T by 

(A(l)C( 2 )et( 3 )) =: -D(l, T) g( 2 , 2 ') g( 3 ', 3 )f ( 2 ', 3 '; 1 '), (50) 

where 

(A(1)A(2)}=:L>(1,2) (51) 

is the static (equal-time) gluon propagator. With the dehnitions of the Grassmann propagator [Eq. (48)] and the 
quark-gluon vertex [Eq. (50)] we can cast the CRDSE for the quark propagator Eq. (49) into the form 

g-i(l, 2 ) = go-i(l, 2 ) + 7 ( 1 , 2 ) - fo(l, 3; 4)g(3, 3')D(4, 4')f (3', 2 ; 4') (52) 

with Qq^ dehned by Eq. (24). Equation (52) is diagrammatically represented in Fig. 1. Alternatively we could have 
put / = ^ in Eq. (47b); this would result in the same CRDSE for the quark propagator as Eq. (52) except that bare 
and full quark-gluon vertex would be interchanged. 

Note that the quark two-point function Q dehned in Eq. (48) is not the true equal-time quark propagator, given in 
the Hamiltonian approach by 

5(1,2) :=i([V;(l),^l(2)]). (53) 

Here the commutator arises from the equal-time limit of the time ordering in the time-dependent Green function. 
Using Eqs. (34) and (35) with f[A] = 1 we hnd from Eq. (53) the relation 

5(1,2) = (e(l)e^( 2 )) + i(A_(l, 2 ) - A+(l, 2 )) = g(l, 2 ) - igo(l, 2 ), (54) 

where we have used the dehnition Eq. (24) of go- For simplicity of notation we will refer to both 5(1, 2) and g(l, 2) 
as quark propagator, since they differ only by a kinematic term go/2 [Eq. (24)]. To exhibit the difference between 
5(1,2) and g(l,2) we consider free Dirac fermions (7 = Tq = 0), for which from Eq. (52) follows g(l,2) = go(l,2) 
and thus from Eq. (54) 

5(1, 2)uee quarks — go(i,2), 

in agreement with Eq. (25). 


^ We denote quark-gluon vertex functions by F, and ghost-gluon vertex functions by F, see Eq. (56) below. 
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-TJSDSDODOB' 


Figure 2. Diagrammatic representation of the CRDSE (60) for the gluon propagator. The small empty dot denotes the bare 
ghost-gluon vertex Fq [Eq. (40)], and the thick dashed line stands for the full (dressed) ghost propagator. The empty square box 
connected to two (amputated) wavy lines represents the variational kernel 2cj, which equals the inverse bare gluon propagator. 



D. Gluon Propagator CRDSE 


The CRDSE for the static (equal-time) gluon propagator Eq. (51) can be obtained by setting f = A in Eq. (47c), 
leading to 

2w(l, 3) D(3,2) - fo(4,3; 1)(^(3) ^t(4) A{2)) - ro(l; 4, 3 )(Ga(3, 4) A{2)) = t(l, 2), (55) 

where we have used Eq. (26). As usual (see Ref. [18]) the full ghost-gluon vertex E is defined by [cf. also the analogous 
Eq. (50) for the quark-gluon vertex] 

{A{ 1 )Ga{ 2 , 3)) = -r( 2 ', 3'; T) D(l, T) G( 2 , 2 ') G(3', 3), (56) 

where G = (Ga) is the ghost propagator. Following Refs. [14, 18] we introduce the ghost loop y by 

2x(l, 2) = fo(3,4; 1) G(4,4') f (4', 3'; 2)G(3', 3). (57) 

In an analogous way we define the quark loop by 

2a(l, 2) = fo(3,4; 1) g(4, 4') f (4', 3'; 2) Q(3', 3). (58) 

Introducing furthermore the gluon energy by 

D(l,2) = iD-i(l,2), (59) 

the CRDSE (55) for the gluon propagator can be cast into the form 

0(l,2)=a;(l,2) + x(l,2) + a(l,2), (60) 

which is diagrammatically represented in Fig. 2. 


E. Quark-Gluon Vertex CRDSE 

The CRDSEs for the quark and gluon propagator, see Figs. 1 and 2, contain the full quark-gluon vertex f defined 
by Eq. (50) and denoted diagrammatically by a full dot connecting a gluon and two quark lines. In the DSEs 
of the functional integral approach in Landau gauge a substantial dressing of the quark-gluon vertex is required in 
order to obtain a sufficient amount of spontaneous breaking of chiral symmetry.® Although in the present Hamiltonian 
approach in Coulomb gauge spontaneous breaking of chiral symmetry is triggered already by the non-abelian Coulomb 
interaction [8, 32, 33] [see Eq. (73) below] the obtained quark condensate, the corresponding order parameter, is far too 
small [32]. It increases substantially when the quark-gluon coupling is included [33]. However, the quark condensate 
obtained in Refs. [33] using a bare quark-gluon vertex is still somewhat too small. At the moment it is unclear 
whether the missing strength of chiral symmetry breaking is due to the use of a bare quark-gluon vertex or due to the 
approximation for the propagators.^ In any case it seems worthwhile to investigate the dressing of the quark-gluon 
vertex; this is given by a CRDSE, which we will derive below. 


® To date, a complete solution of the DSE for the quark-gluon vertex has not been obtained, but models for a dressed quark-gluon vertex 
phenomenologically constructed in accord with the Slavnov-Taylor identity exist, see e.g. Refs. [30, 31]. 

^ In Ref. [33] the variational approach to QCD was formulated in the usual second quantization operator formalism, avoiding the intro- 
duction of Grassmann helds. In that formulation it is convenient to take hrst the fermionic expectation value, which leaves one with 
functionals over the gauge Helds. In the subsequent gluonic expectation valnes in Ref. [33] denominators were replaced by their (gluonic) 
expectation value. This approximation simplihes the analytical calculation but is unnecessary in the present CRDSE approach. 



Figure 3. Vacuum expectation value of two gauge and two Grassmann fields. 



Figure 4. Diagrammatic representation of the DSE (62) for the quark-gluon vertex. 


Putting / = in Eq. (47a) and using (A) = 0 we have 


[Qo-i(l, 4) + 7(1,4)] (C(4) Ct(2) A(3)) + fo(l, 4; 5)(e(4) 7l(5) (2) A(3)) = 0. 


(61) 


The four-point function can be expressed in terms of propagators and vertex functions in the standard fashion by 
Legendre transforming the generating functional of connected Green’s functions to the effective action [18] 

Taking appropriate derivatives of the effective action one finds for the two-quark-two-gluon expectation value 


(A(l)A(2)C(3)?t(4)) = Zl(l,2)Q(3,4) 

+ D{1', 1) D{2', 2) g(3, 3') Q(4', 4)|-f5,^^(3', 4'; T, 2') 
+ r(T,2',5)i:)(5,5')f(3',4';5') 


f (3', 5; T)g(5, 5')f (5', 4'; 2 ') + f (3', 5; 2')g(5, 5')f (5', 4'; T) 


where the two-quark-two-gluon proper vertex is defined by 


1 ) 2 ) 


(5C(4) 6Ai2)SAil) 


A=e=i=o 


This equation is diagrammatically represented in Fig. 3. Inserting this into Eq. (61) the CRDSE for the quark-gluon 
vertex becomes 


f(l,2;3) = fo(l,2;3) 

+ fo(l, 4; 6')g(4, 4')f (4', 5; 3)g(5, 5')f (5', 2; 6)D{6, 6') 

+ fo(l, 4; &)Di4, 4')r(4', 5,3)7^(5, 5')f (6,2; 5')g(6', 6) 

- fo(l, 5; 4)71(4, 4')g(5, 5')f5,AA(5', 2; 4', 3). (62) 

Equation (62) is diagrammatically represented in Fig. 4. Let us also mention that another CRDSE for the quark-gluon 
vertex can be obtained by putting / = AA in the gluonic CRDSE (47c). This results in 

f(l,2;3) = fo(l,2;3) 

+ f (1,4; 6')g(4, 4')f o(4', 5; 3)g(5, 5')f (5', 2; 6)D(6,6') 

+ f(1, 4; 2, 5)g(5, 5')g(4', 4)f o(5', 4'; 3) 

- r,-,ge(l, 2; 4,5)G(5,5')G(4', 4)ro(5', 4'; 3), (63) 


which is diagrammatically shown in Fig. 5. Its new elements are a two-ghost-two-fermion vertex Tqqcc and a four- 
fermion vertex Pqqgg. Both equations are equivalent as long as no approximations are introduced. 
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Figure 5. Alternative CRDSE (63) for the quark-gluon vertex resulting from the gluonic CRDSE (47c). 


-1 


-1 



Figure 6. CRDSE (64) for the ghost propagator. 


F. Ghost CRDSE 

For the sake of completeness we quote here also the CRDSE for the ghost propagator, which was already derived 
in Refs. [14, 18]. To obtain this equation we do not need to explicitly introduce ghost helds. Rather, this equation 
can be obtained already from the operator identity 


Ga{1, 2) = Go(l, 2) - Go(l, 4)A(3) ro(4,5; 3) Ga(4, 2), 

which follows from the dehnition of the Faddeev-Popov operator [Eq. (31)] when this operator is inverted. Taking 
the VEV of this identity and using the dehnition Eq. (56) of the full ghost-gluon vertex we obtain the ghost CRDSE 

G-i(l, 2) = Go-i(l, 2) - ro(l, 4; 3)D(3, 3')G(4,4')r(4', 2; 3'), (64) 

which is diagrammatically represented in Fig. 6. The CRDSEs for the higher n-point functions of the ghost held can 
be derived by representing the Faddeev-Popov determinant as a functional integral over the ghost helds and employing 
the standard DSE techniques, which we are using in the present paper for the quarks and the gluon helds. 

The CRDSEs for the gluon and ghost propagators contain also the full ghost-gluon vertex, see Figs. 2 and 6. The 
CRDSE for this vertex was derived in Ref. [18] and studied in Ref. [34]. It was shown there that its dressing is 
negligible. 

A hnal comment is in order concerning the CRDSEs derived in this section: For purely Gaussian wave functionals 
describing independent quasi-particles these CRDSEs become trivial and are not really necessary, since the higher- 
order correlation functions can be entirely expressed in terms of the two-point functions (propagators) by means of 
Wick’s theorem. However, for interacting theories treated beyond the mean-held approximation non-Gaussian wave 
functionals like our vacuum State [Eqs. (27) and (42)] necessarily emerge. Then the GRDSEs derived above allow us 
to express the various propagators of QGD in terms of the (so far unknown) variational kerneis entering our ansatz 
for the vacuum wave functional. In Sect. IV we will use these equations to express the vacuum expectation value 
of the QGD Hamiltonian in terms of the variational kerneis. In this way these CRDSEs enable us to carry out the 
variational principle for non-Gaussian wave functionals, which are required for interacting helds. 


IV. THE QCD VACUUM ENERGY DENSITY 


With the GRDSEs at hånd we are now in a position to express the vacuum expectation value of the QCD Hamiltonian 
in terms of the variational kerneis. For this purpose we will hrst separate the various powers of the helds in the QCD 
Hamiltonian, so that their vacuum expectation value results in the various n-point functions. 

The QCD Hamiltonian in Coulomb gauge is given by [35] 

i/QCO = -I + i/d»xB“WB“(x) 

+ y’d3crV’”^(x)[-iaia,+/3m]V'”(x)-5 J d^x V'’”^(x)a,A“(x) C„^”(x) 

+ Y j ^^^'^^y^4'V“(x)JA^A'’(x,y) p\y), 


(65) 
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where are the generators of su{N) in the fundamental representation, Ja is the Faddeev-Popov determinant 
[Eq. (30)], and 

is the chromomagnetic held. Furthermore 

FAi^,y) = J <i^zG‘X{x,z){-d^)G'^^{z,y) ( 66 ) 

is the so-called Coulomb kernel, which arises from the resolution of Gauss’s law in Coulomb gauge, and 

p“(x) = /“^M^(x)^^^ + V>’"nx)C„V'”(x) 
is the colour charge density, which we express as 

p(l) = p^(l) + p,(l) =: R{2, 3; 1) A{2) + R{2, 3; 1) V'^(2) ^^(3). 

Here we have introduced the kerneis 

i?(2,3; 1) = - X 2 ) d(xi - X 3 ), (67a) 

R{2, 3; 1) = C 2 m 3 ^3233 ^(xi - X 2 ) (5(xi - X 3 ). (67b) 

For later convenience we rewrite the total QCD Hamilton operator as 

ifQCD = Hym + Hl^ + + ido, Hym =He + Hb + (68) 

and express the various terms in our compact notation. The kinetic (chromoelectric) part of the gauge held then 
reads 

H V-i ^ -T 

^ 2 ^ M(l) ^^24(1)’ 

while the chromomagnetic part can be expressed as 

HB = lj d3x(i?“(x))" =-i24(l)A(l, 2)24(2)+ 1 t(1, 2,3)24(1)24(2)24(3)+ ^T(1,2,3,4)24(1)24(2)24(3)24(4). (69) 

Here we have dehned 


A(l,2) = <5“^“M,,,,a2/(xi-x2). (70) 

Furthermore, the tensor structures T(l, 2, 3) and T(l, 2,3,4) are irrelevant for the present work and can be found in 
Ref. [18]. The non-abelian Coulomb interaction contains a pure Yang-Mills part 


Hé^ = ^JX^PA{l)JAFA{l,2)pAi2), 


(71) 


a fermion-gluon interaction 


= y [JX^PAil)JAFAil, 2) P,(2) + p,{1)Fa{1, 2) pa{2)] , 
and a fermionic Coulomb interaction 


=yP9(i)^^(i.2)p,(2). 


(72) 


(73) 


Finally, the one-particle quark Hamiltonian can be written as 

77d = V'^(l) [ho(l, 2) - J(l, 2; 3) 24(3)] ^(2) = 


(74) 
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Figure 7. Diagrammatic represent at ion of the last term in the r.h.s. of Eq. (77). 


where /iq is defined by Eq. (11) and we have introduced the bare quark-gluon vertex of the QCD Hamiltonian 

J(l,2;3) = gCim 2 (a* 3 )sis 2 < 5 (xi -X 2 ) 5 (xi -X 3 ). (75) 

To carry out the non-perturbative variational approach we now evaluate the expectation valne of the QCD Hamilto- 
nian Eq. (68) in the State defined by Eqs. (27), (41), and (42). We carry out this evaluation up to two-loop level, so 
that the corresponding equations of motion following from the variation of the energy will contain at most one loop. 
Let us emphasize, however, that loops are here defined in terms of the non-perturbative propagators and vertices. 

The magnetic term Hg [Eq. (69)] is insensitive to the quark part of the vacuum wave functional and hence yields 
the same contribution as in the pure Yang-Mills theory. Furthermore, the quark contribution to {Hq^) [Eq. (71)] 
contains more than two loops, which we do not include here. Also the fermion-gluon Coulomb interaction 
[Eq. (72)] yields contributions only beyond two loops and is hence discarded. 

When a Gaussian functional is used for the Yang-Mills sector [see Eqs. (27) and (28)] the cubic term of the magnetic 
energy Eq. (69) does not contribute. Furthermore, the quartic term gives rise to a tadpole in the gluonic gap equation 
[14], which can be absorbed into a renormalization constant. Note also that the contribution of the quartic term to 
the energy vanishes in dimensional regularization. Therefore in the following we will omit the cubic and quartic term 
of the magnetic energy Eq. (69), which then reduces to 

(Hg) = (76) 
where D{2, 1) is the gluon propagator Eq. (51). 

Due to overall translational invariance, the vacuum expectation value of the various terms of the Hamiltonian 
Eq. (65) always contains a diverging factor (27r)^5(p = 0), which is nothing but the spatial volume V. This factor 
disappears when the energy density 

e = {H)/V 


is considered. 


A. Single-Particle Hamiltonian 


The vacuum expectation value of the single-particle quark Hamiltonian Eq. (74) is easily evaluated by means of 
Eq. (35) 


(Hu) = /io(l,2)(V'^(l)^(2)) - J(l,2;3)(^t(i)^(2) A(3)) 

= /io(l,2)[(et(l)^(2)) + A+(2,l)] - J(l,2;3)[(et(l)e(2) A(3)) + A+(2,1)(A(3))]. 

Using (A) = 0 and the definition Eq. (50) of the full quark-gluon vertex T the above expression can be cast into the 
form 


{Hu) = -hoil, 2)[Q{2, 1) - A+(2,1)] - J(l, 2; 3)D(3,3')g(2, 2')f (2', T; 3')Q(1', 1). (77) 

Here the last term arises from the direct coupling of the quarks to the gluons through the bare vertex J(l,2;3) 
[Eq. (75)] in the QCD Hamiltonian. This term is diagrammatically represented in Fig. 7. Besides the bare vertex J 
it contains also the full quark-gluon vertex f, whose CRDSEs have been derived in Sect. HIE [Eqs. (62) and (63)]. 
These relate the full quark-gluon vertex via Eq. (46) to the fermionic variational kerneis. 
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B. Fermion-Fermion Coulomb Interaction 


For the spontaneous breaking of chiral symmetry (SByS) the quark part of the Coulomb interaction Hq [Eq. (73)] 
seems to be crucial. This term alone triggers already SByS [32], albeit not of sufficient strength. On the other hånd 
it was shown in Refs. [33] that the quark-gluon coupling (in Hq,) alone does not provide SB^S, at least within the 
approximation used in Refs. [33].® However, the quark-gluon coupling in Hq, substantially increases the amount of 
chiral symmetry breaking once H'q is included [33]. 

The expectation value of the quark Coulomb interaction Eq. (73) 

E"}? = {H^^) = ^R(3,4; l)i?(5, 6 ; 2 )(^t( 3 ) ^( 4 ) 2 ) ^>^5) V'( 6 )) 

is taken by means of Eq. (36), which yields 

E« = ^R(3,4; l)i?(5, 6 ; 2){ (^+( 3 ) ^( 4 ) 2 ) ^t( 5 ) ^(g)) ( 78 ) 

+ (e^(3)e(4)E^(l,2))A+(6,5) + (C(4)e^(5)E^(l, 2))A+(6,3) 

+ (e^(5)e(6)F^(l,2))A+(4,3) + (et(3)^(6)F^(l, 2))A_(4,5) 

+ [A+(4,3)A+(6,5) + A+(6,3)A_(4,5)] (E.4(l, 2))}. 


Up to two loops in the energy we can replace the Coulomb kernel Fa [Eq. (66)] by its (gluonic) vacuum expectation 
value F = (Fa)- Furthermore, since the Dirac projectors A± are the unit matrix in colour space their contraction 
with the kerneis R [Eq. (67b)] of the quark colour charge density results in the trace of the generators of the gauge 
group, which vanishes.® For this reason the second, fourth, and sixth term in the brackets on the r.h.s. of Eq. (78) 
vanish and we are left with 


E«« ci ^É(3,4; l)i?(5,6; 2)E(1,2) [(^(4) ^t(3) ^(6) C^(5)) + Q(4,5) A+(6,3) - Q(6,3) A_(4,5) + A+(6,3) A_(4,5) 

(79) 

Einally, up to two-loop order in the energy it is sufficient to take the lowest order contribution to the fermion four-point 
function 


(^(4) ^^(3) ^(6) C^(5)) = Q(4,3) Q{6, 5) — (5(4,5) (5(6,3) -I- connected terms. (80) 

Since the quark propagator (5(1,2) is colour diagonal, when Eq. (80) is inserted into Eq. (79) the first term on the 
right-hand side of Eq. (80) gives also rise to a trace over the colour generators and thus to a vanishing contribution 
to the quark Coulomb energy Eq. (79), which then becomes 

E« ^ - y^(3, 4; 1)R(5,6; 2)E(1,2) [Q(4,5) + A_ (4,5)] [Q(6,3) - A+(6,3)]. (81) 


C. The Chromoelectric Energy 


Contrary to the magnetic energy (Hb) and the gluonic Coulomb energy the chromoelectric energy (He) 

does receive additional contributions from the quark sector at the considered two-loop order due to the quark-gluon 
coupling Eq. (43) in the fermionic wave functional Eq. (42). With the explicit form of the vacuum wave functional 
we find from Eq. (29) after an integration by parts with respect to the gluon held 


(He) 


- / VAVS}V£, e"^ 


M(l) 


:. — S*f — \SA 


Ja 


M(i) 


^-Ss-^Sa 


l SS} 1 SSa 1 

r (55/ 1 SSa 1 

yA(l) 2SA{1) 

_,5A(1) 2SA{1)_ 


1/ 5Sa 5Sa\ 1/ SSa S{Sf + S})\ 1 / SS} 6Sf \ 

8\M(1) ,5A(1)/4\M(1) M(1) / 2\M(1) <5A(1)/' 


(82) 


® In Refs. [33] no dressing of the (variational) quark-gluon vertex as described by the CRDSE (62) was included. It is entirely possible 
that when the full dressing of the quark-gluon vertex is included SBxS does take place without including the Coulomb term . In 
faet the results of recent lattice investigations [36] could be interpreted in favour of such a scenario [37]. 

^ This would not be the case within an abelian theory. The arising singular terms can nevertheless be eliminated by an appropriate 
redefinition of the charge operator ^ ^ [ip '^, . 
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To work out the remaining vacuum expectation values we use here the CRDSEs derived in Sec. IIIA. For this purpose, 
using S = Sa + Sf + Sj + [see Eq. (38)] and Sfi/SA = 0, we rewrite the general CRDSE (39) as 


SSa 

SA{1) 



S{Sf + S}) 
SA{1) 



+ ro(3,2;I)(GA(2,3)/). 


(83) 


For the first two terms in Eq. (82) we use this CRDSE with / = 6Sa/SA and obtain 


(55^ \ / S{Sf + S})\ 

,5A(1),5A(1)/ \,5A(1) S Ail) / 


I S^Sa \ 
\M(1)M(1)/ 


ro(3,2;l)^G^(2,3)^^ 


<5(5/ + S)) 8Sa \ 
,5A(1) <5A(I)/’ 


For the last term we can again use Eq. (83) putting / = <5(5/ + 5p/(5A, yielding 


(84) 


biSf-rS)) 8 Sa \ 
<5A(I) M(l) / 


b\Ss + S)) 
M(1)(5A(1) 


+ ro(3,2;I) 


{GAi2,5) 


SiSf + S})\ 

M(l) / 


6iSf + Sf)6iSf + S})\ 

M(l) M(l) /■ 


Inserting this expression into Eq. (84) we obtain 


/ <55.4 <55.4 \ / <55.4 ^(5/ + 5;)\ 

\<5A(I)M(1)/+ \<5A(I) <5A(I) / 

/ ,525 

“ \<5A(I)M(I) 


ro(3,2;I)(|GA(2,3)^^^ 


5(5/+ 5}) <5(5/ +5}) \ 
5A(I) 5A(I) /■ 


Using the above derived expressions we can finally write the Yang-Mills chromoelectric energy Eq. (82) as 


(He) 


1 / å'^S \ 
8\5A(1)5A(1)/ 


1 /5(5/-5p5(5/-5;) 

8\ 5A(1) 5A(1) 


+ lr.(3,2;l)((;,(2,3)+A). 


(85) 


Equation (85) is, so far, exact. Restricting ourselves to the Gaussian ansatz for the gluonic part of vacuum wave 
functional [see Eq. (41)] the various terms can be explicitly calculated. In the last term the definition of the ghost-gluon 
vertex [Eq. (56)] has to be used. One finds then for the chromoelectric energy 

{He) = eI^ + E% ( 86 ) 


where 


= i [D(1,2) - x(1,2)] 0-1(2,3) [0(3,1) - x(3,1)] (87) 

is the contribution arising from the pure Yang-Mills sector, and 

E;Q = -i<T(I,l) + i<T_(I,l) (88) 

is the explicit contribution of the quarks to the kinetic energy of the gluons. In Eq. (88) <7_ is defined analogously to 
<T [Eq. (58)], however, with the bare and full quark gluon vertex, Eq and F, both replaced by 

f _(I, 2; 3) = A+(l, 1') Kil', 2'; 3) A_(2', 2) - A_(1,1') K^il', 2; 3) A+(2', 2). (89) 

[This follows from the terms 5/ — 5^ in Eq. (85).] Using the properties Eq. (13) of the projectors A± and Eq. (24), 
the quantity f _ [Eq. (89)] can be written as 

f _(I, 2; 3) = Qoil, I') fo(l', 2; 3) = -fo(I, 2'; 3) Qo(2', 2). 


(90) 
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D. The Total Energy 

For carrying out the variation of the energy let us summarize the various energy contributions. To the order of 
approximation considered in the present work (two loops in the energy) the total energy is given by [cf. Eq. (68)] 

£^qcd = (Hqcb) — (^^ym) + {Hu) + {Hq), (91) 

with 

{Hym) = {He) + {Hb) + {Hq^). 

Here {Hb), {Hu), {Hq), and {He) are given by Eqs. (76), (77), (81), and (86), respectively. Furthermore, the 
expression for {Hq^) was given in Ref. [14]. As shown in Ref. [25], on a quantitative level this quantity is completely 
irrelevant and will hence be ignored in the following. For subsequent considerations it will be convenient to rewrite 
the energy Eq. (91) in the form 

Eqcd = Eym + Eq , 

where 

Eym = + {Hb) 

is the energy of the Yang-Mills sector, for which we get from Eqs. (76) and (87) 

Eym = i 2) - x(l, 2)] D{2, 3) 1) - x(3,1)] - i A(l, 2)D{2, 1). 

The energy of the quarks is given by 

Eq = {Hu) + E« + E^ (92) 

where the Dirac energy {Hu) = E^'^ + E^'^ was given in Eq. (77). Furthermore, the quark energy Eq. (92) includes 
the quark contribution to the chromoelectric energy [Eq. (88)] as well as the non-abelian Coulomb interaction of 
the quarks [Eq. (81)]. 

The expressions given above for the quark energy can be substantially simplified by noticing that they consist of 
linear chains of fermionic matrices which are connected by ordinary matrix multiplication. Therefore, without loss of 
information we can skip the fermionic indices (but keep the bosonic ones) and assume ordinary matrix multiplication 
for the fermionic objects. This we will do in the rest of this paper’s body. The quark energy Eq. (92) is then given by 

Eq= -Tr[ho(g-igo)] -i?(l,2)Tr[j(l)gf(2)g] 

- iTr[fo(i)gf(i)g] - iTr[gofo(i)gfo(i)gog] 

- ^ E(l, 2) Tr[R(l)(g - igo)i?(2)(g - igo) - lR{l)R{2)] , (93) 

the trace being over fermionic indices only. 

Above we have succeeded to express the vacuum expectation value of the QCD Hamiltonian in terms of the 
variational kerneis (denoted graphically by open square boxes) and the various n-point functions (denoted graphically 
by full dots), the latter being themselves functionals of the variational kerneis through the CRDSEs. In addition, the 
energy contains the bare vertices of the QCD Hamiltonian, denoted graphically by open circles. 

We are now in a position to carry out the variation of the energy. This will result in a set of gap equations, which 
have to be solved together with the CRDSEs. 

V. THE VARIATIONAL PRINCIPLE 

Our trial wave functional [see Eqs. (27), (4I)-(43)] contains three variational kerneis: uj of the Yang-Mills wave 
functional, and Kq and K of the quark wave functional. In carrying out the variations with respect to these kerneis 
we will ignore implicit dependences which will give rise to higher-order loops in the resulting gap equations. This 
implies in particular that we will ignore the dependence of the ghost propagator (and hence of x) on the gluon kernel. 
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Figure 8. Diagrammatic representation of Eq. (94). Small open circles represent the 


“vertices” occurring in the Hamiltonian. 
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Figure 9. Diagrammatic representation of the quark equation resulting from Eqs. (95) and (96). The double line represents 
the Coulomb propagator F = (Fa) [Eq. (66)]. 


as we did already previously in the treatment of the Yang-Mills sector [14]. In the same spirit we will ignore the 
dependence of the gluon propagator on the fermionic kerneis Kq and K as well as the implicit dependence of the 
quark propagator on w with the exception of the free single-particle energy [Eq. (77)], where we will include 
the dependence of Q on the gluon propagator, since this contributes only a one-loop term to the gap equation. The 
explicit derivation of the gap equations is given in App. B. 

The variational equation with respect to w can be combined with the CRDSE (60) as explained in App. B, resulting 
in 


112(1,2) = -A(l,2) + x2(i,2)-2Tr[f(l)g J(2)g] - 2 Tr[f (1) g hogfo(2)g]. (94) 

Here the trace is over the fermionic indices only, as it should be clear from the context. Furthermore, A(l,2) is 
the Laplacian Eq. (70), y is the ghost loop [Eq. (57)], Q is the full quark propagator [Eq. (48)], J [Eq. (75)] is the 
quark-gluon coupling of the QCD Hamiltonian, and Tq [Eq. (46)] is the quark-gluon variational kernel of our trial 
wave functional Eq. (27). Finally, T is the corresponding dressed quark-gluon vertex Eq. (50), which is related to 
the bare one Tq by the CRDSE (62) [or Eq. (63)]. Equation (94) generalizes the gluonic gap equation obtained in 
Refs. [14, 18] to full QCD^° and is diagrammatically represented in Fig. 8. Note that the quarks contribute threefold 
to this equation: i) through the last but one term, which is a quark loop arising from the quark-gluon coupling in 
the QCD Hamiltonian, ii) through the last term, which is a quark loop arising from the free quark energy due to the 
dependence of the quark propagator on the gluon propagator [see Eq. (Bl) below], and iii) through the quark loop cr 
[Eq. (58)] entering the gluon CRDSE (60) for fl. The latter arises entirely from the quark-gluon coupling Tq in the 
QCD wave functional. In App. C the gluon gap equation (94) is used to simplify the expression for the stationary 
energy of the QCD vacuum, which will be needed in future investigations. 

The variation with respect to the biquark kernel Kq leads to the conditions 

K_QhQK+ = 0, A+g/igA_ = 0, (95) 

where h = —5 Eqcy)I5Q is an effective single-quark Hamiltonian 

h = hQ + D{i, 2) [j(i)gf (2) + f (i)g j(2)] 

+ i{fo(i)gf ( 1 ) + f (i)gfo(i)} + i goro(i)gro(i)go 
+ 5 'F(l, 2 )R(l)[g-igo]R( 2 ). (96) 

The quark gap equations are shown diagrammatically in Fig. 9. In the effective single-particle Hamiltonian, ho is the 
Dirac Hamiltonian of free fermions while the remaining terms on the r.h.s. have all the same structure: the quark 
propagator Q (or its free counterpart Qo) is sandwiched by quark-gluon couplings: J is the quark-gluon coupling 
in the QCD Dirac Hamiltonian, Tq and T are, respectively, bare and dressed quark-gluon kerneis of our QCD wave 
functional, and R [Eq. (67b)] is the coupling vertex of the quarks to the Coulomb kernel. 


In Refs. [14, 18] the gluonic Coulomb term [Eq. (71)] was also included, which results in additional terms in the gap equation. 
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Figure 10. Diagrammatic representation of the gap equation (97) for the vector kernel. 


The variation of the energy with respect to the vector kernel K or is carried out in App. B 3. Thereby the 
CRDSE (62) is used to find the variation of the full (dressed) quark-gluon vertex with respect to the kerneis K, ATb 
For simplicity, we quote here the resulting variational equations for K, only in the bare-vertex approximation 
[Eq. (BIO)] 


A± Q{A±fo(l) + 2D{1, 2) [J(2) + hoQfo(2) + fo(2)Qho] }QA^ = 0, (97) 

which is represented diagrammatically in Fig. 10. The kernel K enters here through the bare quark-gluon vertex Eq 
[E q. (46)]. In faet, Eq. (97) can be explicitly solved for K [38]. 

Equations (94), (96), and (97) provide the gap equations for the variational kerneis of our trial ansatz given by 
Eqs. (27), (41) and (42) for the QCD vaeuum wave functional. These equations have to be solved together with the 
CRDSEs for the various propagators and vertices oceurring in the variational equations. In a first step one will use 
the bare vertex approximation, which results in an explicit expression for the vector kernel A'(l, 2; 3) in terms of the 
quark and gluon propagators. Furthermore, one will do a quenched calculation using the gluon propagator (known 
from the present approach to the Yang-Mills sector and also from the lattice calculation) as input for the quark sector. 
Such calculations are presently carried out. 


VI. CONCLUSIONS 

The variational approach to the Hamiltonian formulation of interacting quantum held theories proposed in Ref. [18] 
and developed there for Yang-Mills theory was extended to full QCD. The main feature of this approach is the use 
of CRDSEs to express the vaeuum expectation valnes of powers of held operators (i.e., n-point funetions) in terms of 
the variational kerneis oceurring in the exponent of the vaeuum wave functional. In this way the variational approach 
can be carried out for non-Gaussian wave functionals, i.e., for interacting quantum held theories. To make use of the 
standard DSE techniques this approach requires the use of the coherent fermion state basis of Fock space, which is 
expressed in terms of Grassmann variables. 

Assuming a vaeuum wave functional which contains the coupling of the quarks to the gluons we have derived 
the necessary CRDSEs. By means of these CRDSEs we have expressed the vaeuum expectation value of the QCD- 
Hamiltonian in Coulomb gauge in terms of the variational kerneis of the wave functionals and carried out the variation 
of the energy, resulting in a set of so-called gap equations. These gap equations have to be solved together with the 
pertinent CRDSEs. 

At first sight it seems that the present variational approach is quite eumbersome and less economic than the 
conventional DSE approach in the functional integral formulation of QCD in Landau gauge [1-3]. There one has to 
solve the standard DSEs where the bare vertices are defined by the classical action of QCD. In the present approach 
we have to solve the CRDSEs, which are structurally similar to (and at least as complicated as) the usual DSEs. 
Moreover, contrary to the usual DSEs the bare vertices in the CRDSEs are not known a priori but are variational 
kerneis, which have to be found by solving the gap equations. So it seems that our variational approach is mueh more 
expensive than the conventional DSE approach to QCD in Landau gauge. However, the infinite tower of DSEs has to 
be truncated for practical reasons and there is usually little control over the quality of the approximation achieved. 
Also in our approach we have to truncate the tower of CRDSEs. However, whatever truncation we use, the variational 
principle (i.e., the gap equations) will provide us with the optimal choice of bare vertices for that truncation. We can 
therefore expect that the “bare” vertices of the CRDSEs, i.e., the variational kerneis, capture some of the physics lost 
by the corresponding truncation of the usual DSEs. In faet our “bare” vertices obtained by solving the gap equations 
are not at all “bare” but resemble more dressed vertices of the usual DSE approach [38]. As an illustrative example 
consider the quark-gluon vertex. In the conventional DSE approach in Landau gauge no chiral symmetry breaking 
is obtained when a bare quark-gluon vertex is used in the quark DSE. In our approach we do get chiral symmetry 
breaking even in the bare-vertex approximation. 

In the future we plan to use the approach developed in this paper for a realistic description of the spontaneous 
breaking of chiral symmetry in the QCD vaeuum. Furthermore, we intend to extend this approach to QCD at finite 
temperature and finite baryon density. 
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The present approach is quite general and in principle can be applied to any interacting quantum held theory as 
well as to interacting many-body systems. 
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Appendix A: Diagrammatics 


Propagators: 


Vertices: 


- = g(i, 2 ) 

[Eq. 

(48)] 

-□- = 7(1,2) 

[Eq. 

(44)] 


[Eq. 

(51)] 

— 2cj(l,2) 

[Eq. 

(41)] 

-^- = hoil,2) 

[Eq. 

(11)] 

- = {Fa{1,2)) 

[Eq. 

(66)] 



= T(l,2;3) 

[Eq. (75)] ^ 

^ = i?(l,2;3) 

[Eq. (67b)] 

= fo(l,2;3) 

[Eq. (46)] ^ 

^ = f(l,2;3) 

[Eq. (50)] 

= ro(l,2;3) 

[Eq. (40)] 

1 = r(l,2;3) 

[Eq. (56)] 


Appendix B: Derivation of the Variational Equations 

Below we explicitly carry out the variation of the QCD vacuum energy density with respect to the variational 
kerneis w, Aq, K of our trial ansatz [see Eqs. (41)-(43)] for the QCD vacuum wave functional. 


1. The Gluon Gap Equation 


If the implicit dependence of the ghost and quark loop, x and u, on the gluon propagator is ignored, from the gluon 
CRDSE (60) it is seen that the variation with respect to w can be traded for the variation with respect to D or, more 
conveniently, with respect to the gluon propagator D [Eq. (59)]. 

From the quark CRDSE (52) we have 


resulting in 


6D{3,4) 


hoi2,l)Q{l,l') 


SQ-Hi,2) 

(5D(3,4) 


(5D(3,4) 


Q(2',2) 


fo(l,5;3)Q(5,6)f(6,2;4) 


-ho{2, 1)Q(1, l')ro(l', 5; 3)Q(5,6)f (6, 2'; 4)g(2', 2). 


(Bl) 


(B2) 


Minimization of Eqcb with respect to I?(l,2) then yields the condition 
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-ii?-i(3,l)i?-i(2,4)]i?(4,5)[ii?-i(5,3)-x(5,3) 


1 

+ - 


-1 


1 


-Z?-1(3,1)-x( 3,1) -i?-i(2,3)-x(2,3) --A(2,1) 


2L2“ L2 ^ ’ ^ 2 

- /io(3,4)Q(4, l')f o(l', 5; l)g(5, 6 )f (6,3'; 2)Q(3', 3) - J(3,4; l)g(4, 5)f (5,6; 2)g(6, 3) = O (B3) 

Expressing the gluon propagator Eq. (59) by the gluon energy fl this equation can be simplified to 

12^(1, 2) = -A(l, 2) +x"(l, 2) -2f (5, 6; l)g(6, 3) J(3,4; 2)g(4, 5) -2f (6, 3'; l)g(3', 3)ho(3, 4)g(4,4')f o(4', 5; 2)g(5, 6). 

If we suppress the fermionic indices we arrive at Eq. (94). 


2. The Quark Gap Equation 


The energy depends on the scalar kernel Kq only through the quark propagator Q. Hence the variation with respect 
to Kq and Kq can be carried out as 


SEqcd SEqcd 5g(4,3) ! SEqcd SEqcd (5g(4,3) i 

6Ko{2,l) ~ SQiA,3) SKo{2,l) ~ ’ ^^+(2,1) " ^g(4,3) 6kI{2,1) ~ 

From the quark CRDSE (52) and Eq. (45) we have 


(B4) 


^g(4,3) 

SKo{2,l) 


-g(4,4')A+(4',2)A_(l,3')g(3',3), 


^g(4,3) 

< 54 ( 2 , 1 ) 


-g(4,4')A_(4',2)A+(l,3')g(3',3), 


where we have included only the explicit Kq dependence, since the implicit Kq dependence of Q in the last term of 
Eq. (52) would lead to two-loop terms in Eq. (B4). Defining 


SEqcb 

SQi2,l) 


-h{l,2) 


(B5) 


the stationary condition Eq. (B4) becomes 

A_(l,3)g(3,3')h(3',4')g(4',4)A+(4,2) = 0, A+(l, 3)g(3, 3')h(3', 4')g(4', 4)A_(4,2) = 0. (B6) 

The quantity h{l, 2) [Eq. (B5)] defines an effective quasi-particle Hamiltonian of the quarks. Restricting ourselves also 
up to including one-loop terms in the quark gap equation (B6) only those terms contribute to h{l, 2) which explicitly 
depend on the quark propagator, i.e., the quark energy Eq. (92) of (93). We find then 

h(l, 2) = ho(l, 2) + J(l, 3; 4)g(3, 3')f (3', 2; 4')i4(4,4') + f (1,3; 4)g(3, 3') J(3', 2; 4')14(4,4') 

+ ^{f o(l, 3; 4)g(3, 3')f (3', 2; 4) + f (1, 3; 4)g(3, 3')fo(3', 2; 4) - (f o ^ f _, f ^ f _)} 

+ g^FiA, 4')R(1,3; 4) [g(3, 3') - igo(3,3')] R(3', 2; 4'), 

where f_ is defined by Eq. (89). Using Eq. (90) we recover Eq. (96). 

Equations (B6) are matrix-valued equations. Since A-|_A_ = 0, the expressions on the l.h.s. of these equations are 
manifestly traceless. The relevant information can be extracted by multiplying these equations with Dirac matrices 
and taking the trace. All considerations given above are valid for massive bare quarks. The quark gap equations (B6) 
simplify for massless bare quarks. For instance, multiplying Eqs. (B6) with j3 and taking the trace, thereby using 
/3A-|_ = which is valid for massless bare quarks, we obtain the two conditions 

Tr(g/igA±/3) = 0 


which can be collected in 


Tr(ghg/5) = 0. 
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3. The Equation of Motion for the Vector Kernel 


Finally we derive the equation of motion for the vector kernel K{1,2;3) [Eq. (43)] of the quark wave functional 
[Eqs. (27), (42)]. The energy depends explicitly on the vector kernel K through the bare and full quark-gluon vertex 
[Fq Eq. (46) and F Eq. (50) respectively], and implicitly through the quark propagator Q. Restricting ourselves to 
one-loop terms in the equation of motion we can neglect this implicit dependence in all energy terms except in the 
free single-particle energy [first term on the r.h.s of Eq. (77)]. From the quark CRDSE (52) we get 


sq-HQJ) 

6Kil,2;3) 


'5ro(6,6';8) 
<577(1, 2; 3) 


Q(6',7')7?(8,8')f(7',7; 8 ') 


fo(6,6';8)Q(6',7')i7(8,8') 


5r(7',7;8') 

^77(1,2; 3) ■ 


The variation of the energy with respect to 77(1, 2; 3) yields therefore the following equation of motion (for K'^) 


0 = ho(4,5)g(5,6) 


<5Fo(6, 6 '; 8 ) ^ g,. <5F(7', 7; 8 ') 


L ^77(1, 2; 3) 


<577(1, 2; 3) 


Q(7,4)7^(8,8') 


+ J(4,5; 6)Q(5, 4)0(6, 9) 


1 

+ - 


'<577(1, 2; 3) 

<5Fo(4,5;6)^^^ o gn , 7=^ //< c '5r(7,8; 6) 

<577(1, 2; 3) ^ ,577(1, 2; 3)J 


Q(8,4) 


+ ^go(4,8)fo(8,9; 7)Qo(9,10)Q(10,4) 


(B7) 


For the bare quark-gluon vertex we find from its definition Eq. (46) 


<5ro(4,5;6) 

577(1, 2; 3) 


A+(4,1)5(6,3)A_(2,5), 


^ro(4,5;6) 

577t(l,2;3) 


A_(4,l)5(6,3)A+(2,5). 


(B 8 ) 


On the other hånd, from the CRDSE (62) for the full quark-gluon vertex F, which in the condensed notation of 
Sec. IV D reads 


f (1) = fo(l) + fo(2)gf(l)Qf (3)77(2, 3) + fo(2)gf (3)77(2, 2')77(3,3')r(2', 3', 1 ) - fo(2)gf ,-,^^(1,3)77(2,3) (B9) 

we find the variation of the full quark-gluon vertex F with respect to the vector kerneis 77 and . In taking the 
variation of this equation with respect to 77, 77^ on the right-hand side we can replace the variation of the full vertices 
by those of the bare ones (5f /5K —>■ 5fo/577). This is because the full vertices occur on the right-hand side of Eq. (B9) 
only inside loops and the variation of their dressings would result in more than one loop. Since the variation if the 
bare vertices are explicitly known [see Eq. (B8)], Eq. (B9) provides us with an explicit expression for the variation 
of the full quark-gluon vertex, which has then to be inserted into Eq. (B7). This completes the derivation of the 
variational equations for 77 and C For illustrative purposes we present here these equations also in the bare-vertex 
approximation, replacing the full vertex F by the bare one Fq. The equation of motion (B7) reduces then to 

0 = A_g{77(I, 2) [fo(2)g/io + hogfo(2) + J(2)] 

+ i [2fo(i) + fo(i)go - gofo(i)] }gA+. (bio) 

Using Eq. (90) we can rewrite the last term as 

2fo(i) + ro(i)go - goro(i) = 2(1 - go)fo(i) = 4A_fo(i) 
where we have used the dehnition (24) of Qq in terms of the projectors A±. 


Appendix C: The stationary energy 


For later application we simplify the expression for the energy at the stationary point. For this purpose we multiply 
Eq. (Bl) with 77(3,4) and use the quark CRDSE (52) to find 


<^g-^(i,2) 

577(3,4) 


77(3,4) 


-fo(l, 5; 3)g(5, 6 )f ( 6 , 2; 4)77(3,4) = Q-\l, 2) - go-i(l, 2) - 7 ( 1 , 2). 
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With this relation we obtain from Eq. (B2) 

4) = ho{2, 1 ){Q( 1 , 2 ) - g(l, 1 ') [Qo-i(l', 2 ') + 7 ( 1 ', 2')] Q( 2 ', 2 )}. (Cl) 

Multiplying now the gap equation (B3) with £>(1,2) and using Eq. (Cl) we can express the sum of the Yang-Mills 
energy Eym and quark-gluon interaction energy as 

Eym + 4 '^ = +Eb+ eI^^ 

= i [12(1,1) - x(l, 1)] - hoil, 2)Q(2,1) + hoi2,2')Q{2', 1') [go-'(l', 1) + 7(1', 1)] Q(l, 2), 

where we have used in the last expression the quark CRDSE (52). Adding here also the free Dirac energy E^^ we 
obtain 

Eym+Ed = i [12(1, l)-x(l, 1)] -ho(2, 1) [2g(l, 2)-A+(l, 2)] +ho(l, l')Q(l', 2') [Qo\2', 2)+j(2', 2)] g(2', 2)]. (C2) 

Note that this expression holds only at the stationary point [i.e., for gluon propagators satisfying the gap equation 
(94)]. Formally, the first term on the r.h.s. of Eq. (C2) is the same as the one obtained in Ref. [25] for the pure 
Yang-Mills sector. However, in the present case 12 is the solution of the gap equation (94), which contains the quark 
loop. Also X [Eq- (57)], being a functional of 12 through the ghost propagator, will, of course, be different. 
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